%OFDM系统仿真分析实现
clear all;
close all;
IFFT_bin_length=2048; %IFFT 变换长度
carrier_count=900;%载波数，小于1/2*IFFT_bin_length
bits_per_symbol=2;%多元调制每符号包含的比特数
symbols_per_carrier=10;%发送的OFDM符号总个数

baseband_out_length =carrier_count * symbols_per_carrier * bits_per_symbol;        %基带信号个数(10240)
carriers=(1:carrier_count)+(floor(IFFT_bin_length/4) - floor(carrier_count/2));    %傅立叶变换单元中原始信号放置的傅立叶变换单元号,总长度为实际使用载波数
conjugate_carriers = IFFT_bin_length - carriers + 2;     % 傅立叶变换单元中与原始信号对应的复数放置的傅立叶变换单元号,总长度为实际使用载波数
baseband_out = round(rand(1,baseband_out_length));       %产生baseband_out_length个二进制信号(10240)
convert_matrix = reshape(baseband_out, bits_per_symbol, length(baseband_out)/bits_per_symbol);  %将生成的二进制符号RESHAPE(5120)
for k =1:(length(baseband_out)/bits_per_symbol)                %将 convert_matrix的每一列中的每一行的数值乘以2^ (行数－1)，然后每一列中的各行相加，行数下往上数。
    modulo_baseband(k) = 0;
    for i= 1:bits_per_symbol
        modulo_baseband(k)=modulo_baseband(k) + convert_matrix(i,k)*2^(bits_per_symbol-i);
    end
end

carrier_matrix = reshape(modulo_baseband, carrier_count, symbols_per_carrier)';        %串并转换，每行代表一个OFDM符号
carrier_matrix=[zeros(1,carrier_count);carrier_matrix];                                 %以下4行语句是对每一列进行差分编码
for i = 2:(symbols_per_carrier + 1)
    carrier_matrix(i,:)=rem(carrier_matrix(i,:)+carrier_matrix(i-1,:),2^bits_per_symbol);%R = rem(X,Y) if Y ~= 0, returns X - n.*Y where n = fix(X./Y),fix(a) 为取整。
end
carrier_matrix = carrier_matrix * ((2*pi)/(2^bits_per_symbol));                         %编码后的码字变换到相位
[X,Y] = pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2)));  %极座标到复数座标的转换
complex_carrier_matrix = complex(X,Y);                                                  %组成复信号矩阵
IFFT_modulation=zeros(symbols_per_carrier+1,IFFT_bin_length);                           %傅立叶变换前矩阵初始化
IFFT_modulation(:,carriers)=complex_carrier_matrix;                                     %在傅立叶变换单元中放置复信号
IFFT_modulation(:,conjugate_carriers) = conj(complex_carrier_matrix);                   %在傅立叶变换单元中放置复信号的共轭信号   %
%画出论文中逆傅立叶变换前信号的幅相值特性曲线
figure (1)
subplot(2,1,1);
plot(0:IFFT_bin_length-1, abs(IFFT_modulation(2,1:IFFT_bin_length)),'k*-')
grid on
axis ([0 IFFT_bin_length -0.5 1.5])
ylabel('幅值(伏)')
xlabel('IFFT变换单元')
legend('信号的幅值特性曲线');
subplot(2,1,2);
plot(0:IFFT_bin_length-1, (180/pi)*angle(IFFT_modulation(2,1:IFFT_bin_length)), 'ko')
axis ([0 IFFT_bin_length -200 +200])
grid on
ylabel('相位(度)')
xlabel('IFFT变换单元')
legend('信号的相位值特性曲线');

time_wave_matrix = ifft(IFFT_modulation');                                              %进行逆傅立叶变换获得OFDM的时域信号

print('./public/image/ofdm/ofdm1.png', '-dpng');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v=2;                                                                       %信道的弥散长度
%v=4;
vv=time_wave_matrix(IFFT_bin_length-v+1:IFFT_bin_length,:);                %按书上原理求取循环前缀
sendi=[vv;time_wave_matrix];                                               %加循环前缀
ofdm_modulation = reshape(sendi, 1, (IFFT_bin_length+v)*(symbols_per_carrier+1));  %并串转换
Tx_data = ofdm_modulation;
isi=[0.407 0.815 0.407];                                       %信道特性
% isi=[0.227 0.460 0.688 0.460 0.227];
Tx_data=conv(Tx_data,isi);                                       %发送信号与信道冲激响应卷积
Tx_data =Tx_data(2:length(Tx_data)-1);
%Tx_data =Tx_data(3:length(Tx_data)-2);

time_wave_matrix = real(time_wave_matrix');                %取实部
%画出无加循环前缀OFDM符号的时域波形曲线
figure (2)
subplot(2,1,1)
plot(0:IFFT_bin_length-1,time_wave_matrix(2,:),'k')
grid on
ylabel('幅值(伏)')
axis ([0 IFFT_bin_length -0.07 0.07])
xlabel('时间(ts)')
legend('一个OFDM符号时域波形')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ofdm_modulation = reshape(time_wave_matrix', 1, IFFT_bin_length*(symbols_per_carrier+1));
%画出无加循环前缀OFDM符号的频谱    直接对时域的相频特性和幅频特性做变换。
temp_time = IFFT_bin_length*(symbols_per_carrier+1);
symbols_per_average =ceil(symbols_per_carrier/5);
avg_temp_time = IFFT_bin_length*symbols_per_average;
averages = floor(temp_time/avg_temp_time);
average_fft(1:avg_temp_time) = 0;
for a = 0:(averages-1)
    subset_ofdm = ofdm_modulation(((a*avg_temp_time)+1):((a+1)*avg_temp_time));
    subset_ofdm_f = abs(fft(subset_ofdm));
    average_fft = average_fft +(subset_ofdm_f/averages);
end
average_fft_log = 20*log10(average_fft);
subplot(2,1,2)
plot((0:(avg_temp_time-1))/avg_temp_time, average_fft_log,'k')
axis([0 0.5 -40 max(average_fft_log)])
ylabel('幅值(dB)')
xlabel('频率（fs/2＝0.5（HZ）)')
legend('OFDM信号频谱')

print('./public/image/ofdm/ofdm2.png', '-dpng');


Tx_signal_power = var(Tx_data);  %求取发送信号的功率为加高斯白噪声作准备
for SNR=1:50
linear_SNR = 10^(SNR/10);
noise= Tx_signal_power/linear_SNR;
noise= sqrt(noise);
noise = randn(1, length(Tx_data))*noise;  %
Rx_Data = Tx_data + noise;           %加高斯白噪声
Rx_Data_matrix=reshape(Rx_Data,IFFT_bin_length+v,symbols_per_carrier+1);   %串并转换
Rx_Data_matrix=Rx_Data_matrix(v+1:(IFFT_bin_length+v),:);                  %去循环前缀
Rx_spectrum = fft(Rx_Data_matrix);                                         %傅立叶变换
Rx_carriers = Rx_spectrum(carriers,:)';                                    %取有用信号
Rx_phase = angle(Rx_carriers)*(180/pi);                                    %计算相位
phase_negative =find(Rx_phase < 0);                                        %找出相位为负的位置
Rx_phase(phase_negative) = rem(Rx_phase(phase_negative)+360,360);          %把相位转换到0－360度
Rx_decoded_phase = diff(Rx_phase);                                         %相位差分
phase_negative = find(Rx_decoded_phase< 0);                                %找出相位为负的位置
Rx_decoded_phase(phase_negative) = rem(Rx_decoded_phase(phase_negative)+360,360);%把相位转换到0－360度
base_phase = 360/2^bits_per_symbol;                                        %求取相位基准
delta_phase = base_phase/2;                                                %相位辨别门限
Rx_decoded_symbols = zeros(size(Rx_decoded_phase,1),size(Rx_decoded_phase,2));
%解码:
for i = 1:(2^bits_per_symbol - 1)
    center_phase = base_phase*i;                           %相位判决中间值
    plus_delta = center_phase+delta_phase;                 %相位判决上限
    minus_delta = center_phase-delta_phase;                %相位判决下限
    decoded = find((Rx_decoded_phase <= plus_delta) & (Rx_decoded_phase > minus_delta));   % 找出所有符号中在各判决区域内的点
    Rx_decoded_symbols(decoded)=i;                         %将以上各点译码为1～2^bits_per_symbol - 1，没变者为初始值0
end
Rx_serial_symbols = reshape(Rx_decoded_symbols',1,size(Rx_decoded_symbols,1)*size(Rx_decoded_symbols,2));   %并串转换
for i = bits_per_symbol: -1: 1                             %按模2取余变为2进制
    if i ~= 1
        Rx_binary_matrix(i,:) = rem(Rx_serial_symbols,2);
        Rx_serial_symbols = floor(Rx_serial_symbols/2);
    else
        Rx_binary_matrix(i,:) = Rx_serial_symbols;
    end
end
%计算误码率
baseband_in = reshape(Rx_binary_matrix,1,size(Rx_binary_matrix,1)*size(Rx_binary_matrix,2));
bit_errors = find(baseband_in ~= baseband_out);
bit_error_count = size(bit_errors,2);
snr(SNR)=bit_error_count/baseband_out_length;
end
%画出误码率曲线
figure (3)
plot(snr);
ylabel('误码率')
xlabel('信噪比（dB）')
legend('误码率曲线')

print('./public/image/ofdm/ofdm3.png', '-dpng');
